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We study two free energy approximations (Bethe and plaquette-CVM) for the Random 
Field Ising Model in two dimensions. We compare results obtained by these two methods in 
single instances of the model on the square grid, showing the difficulties arising in defining 
a robust critical line. We also attempt average case calculations using a replica-symmetric 
ansatz, and compare the results with single instances. Both, Bethe and plaquette-CVM 
approximations present a similar panorama in the phase space, predicting long range order 
at low temperatures and fields. We show that plaquette-CVM is more precise, in the sense 
that predicts a lower critical line (the truth being no line at all). Furthermore, we give 
some insight on the non-trivial structure of the fixed points of different message passing 
algorithms. 
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I. INTRODUCTION 


Many of the thermodynamic properties of the Random Field Ising Model (RFIM) remained 
controversial for more than 40 years. In particular the existence or the absence of a spin glass 


phase in the model attracted a lot of attention 


JJ-rI] . Only recently we found confirmation that t 


Uis 


spin glass phase is possible only if the magnetization of the system is constrained to be fixed 8|]. 
Without this constraint, the RFIM possesses only a paramagnetic-ferromagnetic transition when 


d > 2 or no transition at all in d < 2 


a-lnt] 


These results bring new questions into the field. For example: Why standard perturbation 
theory fails? Could the spin-glass phase, deduced for the model with fixed magnetization in 
the Bethe lattice Q], have some role in the dynamical behavior of models in finite dimensions? 
Do message passing algorithms always converge in lattices of finite dimensions as occurs in the 

121 ? 


ferromagnetic Ising model 


On the other hand, our current knowledge about the model makes it a good starting point to 
improve the comprehension about the applicability and the capabilities of techniques that have 
being used to study this and other disordered systems. Of special interest is the analytical solution 
by Mezard and Parisi 13, Q of the Viana-Bray model [l^) within a Replica Symmetry Breaking 
ansatz. Since then, the field has progressed very fast. First, the ansatz was rapidly extended to 
other models 


16Hl8t|. Then, it was soon recognized the connection between this approach and 
message passing algorithms, the well-known Belief Propagation (BP) algorithm [l^ corresponds 


to the Bethe approximation of the free energy of a specific model 


2C|. 


To go from the Bethe approximation to methods that considered loops in the interaction net¬ 
works turned out to be a more difficult task 


from Yedidia and co-workers 


2ll-|28l|. An important step in that direction came 


29(1 that described how to generalize the Cluster Variational Method 


(CVM) of Kikuchi [30|. The minimization of the CVM free energy can be achieved by the use 


of a Generalized Belief Propagation (GBP) algorithm 
Replica Symmetric (RS). 


29( 1. although the solution found is always 


More recently it was possible to merge the CVM with the Replica Symmetry Breaking (RSB) 


ansatz 


311 1. Although, technically involved, the approach allowed the study of the Edward-Anderson 


(EA) model in two dimensions. In this system it was possible to unveil the connection between the 
phase transition predicted by the average case scenario and the properties of the Generalized Belief 
Propagation algorithm in two dimensional lattices 321 • Running standard BP for the Bethe approx¬ 


imation in EA 2D one finds a paramagnetic solution at high temperature, as expected. However, 
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decreasing temperature BP finds, not one, but many fixed points with non-zero local magnetiza¬ 
tions. Suggesting then, a transition from a paramagnetic to an spin glass phase that although 
is not expected in thermodynamical grounds, is nevertheless predicted within the approximation 
and connected with the performance of message passing algorithms 




33( 1 ■ Furthermore, these 


low temperature fixed poi nts are also correlated with the metastable states of the Monte Carlo 
dynamics of the model 341. 


In this work we present the equations describing the RFIM within the CVM at the RS level. 
The goals are, on one side, to compare the predictions that can be derived from the equations with 
our current knowledge about the model. This should shed light on the limitations and advantages 
of the approximation itself. On the other, to extend the connection already established within the 
EA model between these average case equations and message passing algorithms. In particular, we 
will focus on whether the CVM at the Bethe and plaquette level also predict a spin-glass phase, 
and if this is connected with the behavior of message passing algorithms. 

The rest of the work is organized as follows. In the next section m we present the model. In 
section m we obtain message passing equations for the single instances of the model, both at the 
Bethe (BP) and the plaquette level of CVM (GBP). The phase diagram obtained by implementing 
these belief propagation algorithms is discussed in IIVI In section |V] we derive average case predic¬ 
tions at replica symmetric level for both approximations. We summarize and discuss our findings 
and conclusions in lYD 


II. THE MODEL 

The Random Field Ising Model is a natural extension of the standard Ising ferromagnet to 
consider the disorder of real solids. Due to a number of reasons, spins in a crystalline lattice 
are under the influence of small local magnetic fields; intensity and direction of these fields varies 
rapidly from one site to another uncorrelatedly. This situation can be modeled in a simplified way 
by considering the following Hamiltonian; 

nia) = CTiaj - ^ hiai 

{*J> * 

where the first sum runs over all couples (i,j) of neighboring spins (first neighbors on the lattice) 
and the second over all sites (spins). We will deal in this work with N = L x L spins in a 
bidimensional square lattice, but extensions to more dimensions is conceptually straightforward, 
although probably numerically cumbersome. The magnetic exchange constants between spins is 
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fixed to J > 0 (more precisely J = 1) and the spins ai = ±1 are the dynamic variables. The 
disordered local field hi are a set of random variables generally drawn from a zero-mean bimodal 
or Gaussian distribution. In our case, a bimodal distribution will be used: 

Ph [hi] = \ {5{h, -H) + 6{h, + H)) 


Local fields introduce what is called quenched disorder 


351]: fields and spins are considered to 


vary in completely different time scales in real solids. Thus, an instance of this model corresponds 
to a particular realization of the fields. We will use the field intensity H as one of the model 
parameters. 

The statistical mechanics of the RFIM model, at a temperature T = 1//3, is given by the 
Gibbs-Boltzmann distribution 


P{a) = 


z 


where Z = i 


-mp) 


is a normalization constant, called partition function since it depends on the temperature T = 1//3 
and all other relevant parameters of the model, like H for instance. 

Finding the exact free energy of a particular sample of our model is quite difficult, because the 
partition function involves an exponential number of computations, being only accessible to very 
small systems. However, the interesting limit is the thermodynamic one {N —>■ oo), and from a 
physics perspective we are mostly interested in a description of the ensemble of instances, not in 
a particular one. Using the replica method for models in fully connected or very sparse random 
topologies (see [35| for a primer), the structure of equilibrium states for the average case can be 
described in a well defined hierarchy of approximations. The whole approach is based on the idea 
that the free energy is a self-averaging quantity, so averaging over disorder must match results for 
real samples when N —>■ oo. An exact solution to disordered models in finite dimensions, however, 
remains immune to analytic results, and usually is treated by simulations or approximations, like 
the ones described next. 


III. GBP EQUATIONS FOR PLAQUETTE-CVM 


The fundamental idea behind Generalized Belief Propagation 
Method 




29(, or Gluster Variational 


37], or Kikuchi approximation |30| to the free energy of a model is to replace the 


exact functional free energy F[P(fTi, ct 2 , ..., uat)] by an approximation consisting on the sum of 









5 


local (small regions) free energies: 


F ~ -Fcvm = ^ CRFR[Pji{aR)]. 




Since the set of regions is in principle arbitrary, some regions might contain smaller regions, and 
therefore an appropriately chosen counting number cr is used to avoid overcounting free energy 


contributions 


3fil |. Only in very particular cases a partition like this remains exact (like in t 


ue 


Bethe approximation on tree like topologies) or even an upper bound of the real free energy 


361]. 


Ideally the distributions minimizing Fqym should coincide with the marginals of the exact 
Boltzmann distributions Pr{(7r) = *^ 2 ) ■ ■ ■, c^at), but generally (and hopefully) they 

will only be an approximate of them, and therefore are usually called beliefs and represented by 
bR{aR), leaving PR{crR) for the exact marginals. 



Plaquette 


Link 


Spin 


FIG. 1. Regions for CVM approximation 


Though very general, there are specihc ways to define the set of regions TZ for approximating 
the free energy. According to the CVM prescription [^ . for instance, this set is formed from a 
primary set of larger regions TZq. Then adding to them the set TZi, formed by the intersections of 
every pair of regions in TZq. Then adding the intersections of regions in TZi and so on and so forth 
until no intersections are left, obtaining TZ = TZq U TZi U TZ 2 ■ ■ ■■ 

Bethe approximation is just a particular case of CVM where the larger regions considered are 
the interacting pairs of spins (links). In this work we use Bethe approximation, and, to go beyond, 
we also use the approximation starting from plaquette regions. A plaquette region is formed by the 
four spins and links lying on a basic square cell of the lattice (see Fig. [T]), and has the advantage 
over Bethe approximation of including the shortest loops in the graph. It is known that loops 










6 


are the cause of failure of Bethe approximation. Two adjacent plaquettes overlap in a link region, 
formed by two neighbor spins and the interaction between them. Two link regions with a common 
spin overlap in a spin region. Therefore, our approximation contains contributions of free energies 
from all plaquettes (P), links (L) and sites (i) (Fig. [T|). The free energy approximated with this 
set of regions is usually called Kikuchi approximation [30|, l37| : 


FcvM[{bp},{bL},{bi}] = '^Fp[bp{aa,(Tb,ac,ad)] FL[bL{<Ta,crb)] +'^ Fi[bi{ai)] 

P L i 

Minimization of this free energy with respect to the beliefs {&p}, has to be done 

while keeping the consistency among them, in the sense that distributions on larger regions have 
to marginalize onto the smaller regions that are part of them. This constrained minimization is 
implemented via Lagrange multipliers Mp^p and mp^i enforcing plaquette to link and link to spin 
consistency respectively. The self consistent equations for the Lagrange multipliers are interpreted 
as message passing equations. The message from link (i,j) to spin, say, i is updated according to: 


PePKij)] LeC{j)\{i,j) 

where 'P[{i,j)] is the set of regions of whom (i,j) is a child. In Fig. [23 it is shown the configu¬ 
ration of all messages involved in updating The interactions appear as Uj) = 

exp(—/3Pjj(iTj,(Tj)), where Eij{ai,aj) = — Jaiaj — hjaj. The product over link-to-spin messages 
runs on C{j), which is the set of links having the site j as a child, excluding {i,j). In the case 
shown, C(j) \ {i,j) includes {h,j), {k,j) and (/, j). 

For messages from the square plaquette A = {ijkl) to link {i,j) (see Fig. [20 we have the 
update equation ([2]) that is similar in structure. This time, however, on the left hand side appears 
the product of two link messages that are internal to the plaquette. Normally it is advisable (for 
convergence issues) to update these first before updating j). 




n n Mp^p{ap)\ j] 


yL&Cp[A]\(pj) P&V[L]\A 


)( 


s •>^1 


n 


fklp^nis^n) 


( 2 ) 


ns{fc,d LaCs['n\/\L^Cp[A\ 


Here, Tp[H] is the set of links in plaquette A, Cs{n\ represents the links that are parent to spin n 
and V{L\ is the collection of parent plaquettes to L. 

Messages are the usual representation of the Lagrange multipliers in the context of belief prop¬ 


agation algorithms. However, physically it is nicer to think in terms of cavity fields 


Q,Q, 


acting 
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(a) Link to spin (b) Plaquette to link 

FIG. 2. Schematic representation of message update equations. Green messages are updated using black 
and red ones. Red messages appear on the LHS of the respective equation. 


from one part of the system onto another. In terms of these fields, messages are recovered as 


= exp/3uLicri = exp ,0(f7f7icrj + Uiai + UjOj) 


(3) 


Only one parameter u ta is enough to represent mi^i messages, while for we need three: 


Upij, upi and upj 


31 


M, 


381 ] ■ The update equations for the cavity fields are obtained with little 


effort. Below we write them using the labeling of Fig. [2j 


= —arctanh[tanh (3J tanh/3/ij] + up^^i + upg^i 

P 

where 

hj hj Ui^j + Uj^^j Ufi^j Upj^^j Up^^j 

J J ~\~ UU 

Fields define a new effective hj and coupling constant J. For plaquette-to-link fields we get 


(4) 


lnK(_l,l)i^(l,_l) 

1 , K(1,1)K(1,-1) 


( 5 ) 


4/3 7^(-l,l)K(-l,-l) 


1 1 iF(l,l)iF(-l,l) 

“ 1 “ A 


4/3 i^(l,-l)iF(-l,-l) 
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where 


K{ai,aj) = ^ exp/3 + Uc^kicrkcri + UB^jkfyj^k + + 

(o-fc.o-i) 

The set of self consistent equations defining messages can be solved by a fixed point iteration. 
This procedure is known as a message passing algorithm. Standard Belief Propagation (BP) is 
just a special case for the Bethe approximation. Many details about the implementation of BP 
and GBP can be found in 33|. If the algorithm happens to find a fixed point (convergence is not 


guaranteed), you can recover the local approximated marginals 6 _r(ct/j) in terms of the messages, 
just like with any standard Lagrange multiplier minimization. 

Now that we have presented in broad terms the essential features of BP and GBP, we can move 
to discussing the results when rnnning these fixed point methods on singular instances of 2D RFIM. 


IV. RESULTS FOR SINGLE INSTANCES 

Region graph approximation to the free energy and the corresponding message passing algo¬ 
rithms give numerical (estimate) values to all thermodynamic quantities like free energy, energy 
or entropy. Usually these approximations are compared with Monte Garlo simulation results, or 
exact predictions (when available). We are most interested in the phase diagram of the model. 

Two relevant parameters characterize the random field Ising model: the temperature T = 1/13 
and the external field intensity H. Local fields hi = Hhi are scaled by H, and hi G {—1,1} with 
equal probability. The model with LI = 0 is no other than the usual Ising ferromagnet. 


Phase Diagram 

We know in advance the dominant thermodynamical phase of some particular points in the 
(r, H) plane. For instance, at any finite field intensity H, and high enough temperature T, the free 
energy minimization is dominated by the entropic term, and therefore spins are mostly uncorrelated 
(true when T —)• oo) with vanishing local magnetizations, rrii ~ tanh/3/ij ~ I3H = H/T, and with 
global zero magnetization M = ^ ~ H/T given the random origin of the fields hi. 

At zero temperature, the free energy is dominated by the minimization of the energetic term, 
and at least for H > 4, the interaction with the 4 surrounding spins at any given site cannot 
overcome the interaction with the local external field. In such case, the ground state of the system 
is trivially ViCTj = Sign(/ij). The system is frozen and again global magnetization is zero (except 
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for 1/y/N fluctuations). Let’s call this phase paramagnetic in spite of spins been magnetized. In 
this sense, paramagnetic refers to the fact that no long range order exists, and the frozenness of 
the spin is not a collective phenomena but a result of their interaction with external local fields. 

Therefore, in both extremes, high T and high H, the system is paramagnetic, and this is 
correctly reproduced by both, Bethe and CVM approximations, on single instances (see below). 
To signal the onset of long range correlations we have used the global magnetization as the order 
parameter and a threshold value \M\ > Mt^. The value of this threshold can be fixed according to 
the following reasoning. For a given temperature, long range order, if present, is destroyed for large 
H. In this regime, magnetization is of order m ~ = -^. If the system remains paramagnetic 

when decreasing H, magnetization should decrease too. On the contrary, if for lower H states with 
m S> are found, it can be interpreted as the emergence of order in the system. In practice, a 
value of Mth = ^ fits known results for FI = 0, namely the para-ferro transition point of the Ising 
model. 

To unveil the phase diagram in the thermodynamic limit from simulations on finite size instances 
we considered that, for a given value of the H field, the mean value of the critical temperature 
over many samples depends on the system size, Tc = Tc{L). In order to obtain an estimate of the 
value in the L —)• oo limit, we propose \TJL) — Tc(oo)| oc —^ — and extrapolate (see figure [3]). 

Summarizing, in figures H] and [S] we show the average phase diagram over many samples of 
RFIM with different system sizes and the critical line after extrapolation. Our Bethe and CVM 
calculations predict the appearance of long range correlation in the magnetization of spins at low 
temperatures and low external fields. This phase (thermodynamically incorrect for L —>• oo) will 
be called ferromagnetic. In both cases the para-ferro border starts at the corresponding critical 
point of the ferromagnetic (non disordered, H = 0) model. This transition is at odds with the 
rigorous result proving that in d < 2 even the smallest (random) external field should destroy all 


chances of a long range order 


Q, 11 , 3- 


It is not a surprise that mean field approximations stabilize non relevant phases, as it does in 
the Edwards-Anderson model [3 . 3- This is a general price we pay for implicit factorization at 
certain levels of the correlations when writing the free energy approximated by regions. 

In the CVM plaquette approximation, we faced severe convergence problems that we solved 
partially by fixing the gauge invariance of message passing equations [3^ as explained in appendix 


IVIIl In the following, CVM approximation refers to this gauge fixed implementation of the message 
passing. However, even in this case, there is an island of non convergence for low temperatures 
and intermediate values of field (see B- 
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Extrapolating Tc 
Kikuchi(CVM), H=1.0 



FIG. 3. For a fixed field intensity, different T^. values are found depending on the size of the lattice. Here 
we assume \Tc{L) — Tc(oo)| oc make a linear extrapolation for L —>■ oo. This procedure is repeated 

for each H value in the CVM and Bethe approximation. 



(a) Critical lines for L=16,32,64. (b) Extrapolated Tc{L) for L —>■ oo. 

FIG. 4. Bethe(BP) approximation. Note in both plots that for H ^ 0 the result for the Ising magnet is 
recovered. 

This region of non convergence grows with the system size. In Fig. [5b] we show, for different 
values of L the border of the area where the probability of convergence drops abruptly. As L 
increases, the non convergence island apparently spreads all over the ordered phase. To be more 
quantitative on this point we plotted for a fixed T = 0.5 in Fig. [6a|the value of H = H^c{L) at 
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(a) Critical lines for L=16,32,64,128. (b) Extrapolated L —>■ oo critical line. 

FIG. 5. Kikuchi(CVM) approximation. In part (b), for {T, H) pairs inside the region limited by dashed lines 
the algorithm does not converge. The non convergence region seems to spread over the whole ferromagnetic 
phase as L increases. 

the frontier for each L. The functional dependence of this curve is very well fitted with the law 
H{L) = CL~^. This suggests that indeed, for large enough lattices the non convergence region 
moves down to H = Q. 

Also, for a fixed = 1.5, we fitted in Fig. [6b] the temperature at the convergence border. The 
form of T^c{L) in this case agrees with Tc — CL~^. The limit lim T^q{L) = Tc = 1.28 is close to 

L^OO 

the average critical temperature found by GBP for the same field intensity, Tc{H = 1.5) = 1.40. 
These results enable us to speculate that in the thermodynamic limit the difficulties in convergence 
would cover all the unphysical ordered phase. Moreover, this also suggests that at least for this 
model the convergence of GBP is linked to the para-ferro phase transition and may be useful in 
defining its location. 


Multiple solutions: a one case study 


For many years there was a strong debate about whether a spin glass phase was present in 
an intermediate area of the phase diagram just around the border of the para-ferro transition 


Q, 


lOj, I 111 . [39(1 . Although we now know that long correlation order is impossible in this ca^ such 
belief was supported by perturbative and non perturbative replica field theory resultsjl, 40]. On 
the other hand, although numerical Monte Garlo simulations suffer from a strong slowing down 
near the transition, no evidence of a SG phase is found in either the bimodal or the Gaussian 


version of the RFIM 


4l|-[43|. 
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Fitting non convergence H,sjq(L) 

Kikuchi(CVM), T=0.5 
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(a) For a fixed T, H{L) = ChL-*=^ fits well the field 
value at the lower border of instable regions. 


Fitting non convergence T|vjq(L) 

Kikuchi(CVM), H=1.5 



(b) For a fixed H, T{L) = T^- CtL-^^t fits well the 
temperature value at the right border of instable regions. 


FIG. 6. Large L behavior of the convergence frontier (see Fig. I5bl) . For L —>■ oo, the approximately flat 
lower border of the instable region approaches the line H = 0. On the other hand, the right border seems 
to approach the critical line. 


Solution Multiplicity Bethe 

T=0.5 



Solution Multiplicity Bethe 

T=1.5 



FIG. 7. Number of solutions found for two samples of different size as a function of field intensity. Multiple 
solutions appear for H < He, within the ordered phase. For small H the number of solutions reduces as for 
H = 0 the solution is unique. The number of stable fixed point for the BP equations increases with the size 
of the system and decreases with temperature. 


It is consistent with such expectation the fact that Bethe approximation finds many different 
solutions in the vicinity of the para-ferro transition. They are discerned easily by their different 
magnetizations. In Fig. [7| we show the number of different solutions found in a system of = 
32 X 32 and Af = 64 x 64 at T = 0.5 and T = 1.5. The message passing is started from random 
initial conditions 10^ times for each value of the external held intensity. Changing the external held 
means varying only its intensity while keeping directions on each site. As expected, the larger the 











13 



FIG. 8. These are four snapshots of the magnetization of a 32x32 system. Starting from different initial 
conditions for the same field values, the BP algorithm finds different solutions characterized by domain 
formation. 

system the more solutions are found by the algorithm. Also, for lower temperatures the number of 
solutions increase, as new local energy minima get stabilized. Notice that multiplicity of solutions 
appears for H ~ Hc{T) and below, within the long range order phase. As a consequence, a proper 
dehnition of the critical held/temperature with a threshold in the global magnetization is even 
more difficult than discussed previously. 

Different Bethe solutions are not completely uncorrelated. They differ in blocks of spins that 
switch directions from one solution to the other, as shown in Fig. [8l It has been shown for other 
models, like the Edwards-Anderson in 2D, that these solutions might be connected with the regions 
of the phase space in which the Monte Carlo dynamics spend more time, or say in other words 
where the Boltzmann measure concentrates Q- It is to be checked in future works whether this 
property holds for the RFIM too. 

Furthermore, we wonder if this multiplicity of states is truly a sign of a spin glass like behavior 
for hnite lattice sizes, or also an artifact of the Bethe approximation. More precisely, we want to 
check if it is also a feature of plaquette-CVM approximation. Surprisingly it is not. In most systems 
studied, only one GBP hxed point was found. In some of them, two or at most three solutions were 
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found. When more than one solution is found, they are still related by the flipping of entire blocks 
of magnetizations. The drastic reduction of fixed points for plaquette-CVM is conspicuous, and 
might be interpreted as a sign of the non thermodynamic relevancy of the many Bethe solutions. 
This is in accordance with the general belief that CVM yields more accurate results than Bethe 
approximation. 


Solution Multiplicity Bethe Solution Multiplicity Bethe 

64x641=0.5 64x64 T=1.5 



FIG. 9. For a given 64x64 system, different magnetizations are obtained for each H value when running BP 
with various initial conditions. The solution with the lowest free energy is drawn with a solid line. This line 
should be useful in defining the transition more accurately. The He values correspond to the temperature 
of each plot and are taken from the results shown in Fig. |4l 


This detailed study of the Bethe and CVM approximations shows the difficulties affronted in 
defining a critical line in the random field Ising model. At variance with the Ising ferromagnet, 
where the Bethe and CVM approximations give a clear zero magnetization in the paramagnetic 
phase[l^, here there are non zero magnetizations at every temperature and field, and the onset of 
long correlation is further obscured by the appearance of many solutions. It would be more justified 
if one could define the critical line studying a threshold magnetization value for the solution with 
lowest free energy at every (T,H). The different magnetization values for a 64x64 sample are 
plotted in Fig. [9l Solutions with the lowest free energy are joined with a solid line. Applying 
the threshold criteria to this line should improve the quality of the results. However, this requires 
running message passing thousands of times from random initial conditions for each system, and 
then the average over thousand of systems. This is a highly resource-consuming task. Instead, we 
turned to the prediction of a critical line by average case calculations of both Bethe and plaquette- 
CVM approximations. 
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V. AVERAGE CASE PREDICTIONS 


There are two different procedures to analyze the average case solution of CVM approximations, 
one more formal relying on the replica method 


based on population dynamics 


3ll | and quite involyed, and another more intuitiye 


32l |. In the case of the Bethe approximation both procedures are 


equivalent, the latter being a numerical solution of the integral equation appearing in the former 
13l |. A detailed explanation of both procedures can be found in 3^, and we will only describe 
them shortly. 

Replica-CVM is a formal way of averaging the disorder in the partition function of a region graph 
approximated disordered model, based on the replica trick 311]. As usual in replica method, the 


average over the disorder couples the replicas, and some ansatz is needed to take the n ^ 0 limit. 
Generalizing the approach in [^, 45| to the plaquette-CVM case, this ansatz is a parametrization 
of messages m and M in terms of two types of field distributions q{u) and Q{U, ui, M 2 ) respectively; 


m{ai) (X j duq{u) exp 


af 


a=l 


M{ai,aj)(x / dUduidujQ{U,Ui,Uj)ex.p 


pu Y Y 


CTn 


( 6 ) 

(7) 


a=l a=l a=l 

Minimizing the replicated free energy in terms of these distributions, we obtain two self consistent 
equations for q{u) and (5(?7, mi, M 2 ): 

k p 

q{u) = j Y[dqiY[dQa{5{u-u{H^)))h 

i a 

R(U, Ua, Ub) = J duidujQ(U, Ui,Uj)q{ua - Ui)q{ub - Uj) = 

. K p 

= / Y\dqiY\dQa{d{U -U{#))6{Ua-Ua{#))d{ub-Ub{#)))h 

i a 

with u{#) and U{#) functions defined in (|3|) and ([S]) and 


( 8 ) 


dqi = q{ui)dui 


dQa ^ QaiU,Ui,Uj)dUduiduj 


All the thermodynamics is encoded into the functions q and Q. The solution of this system is 
technically involved. Standard population dynamics is not easily applicable because of the need 
of numerical deconvolution methods (like Fourier transform) to extract Q from the convolution in 
the left side of Eq([5]). In the case of zero external field and high T, however, there is the trivial 
solution q{u) = S{u), and Q(U, m i , M 2 ) = Q{U) 6 {ui) 6 {u 2 ), with Q{U) satisfying a simplified version 
of the integral equation ([8]) (see [3^ for details). 
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We reproduce the stability analysis of the paramagnetic solution for the RFIM, knowing in 
advance, that q{u) = 5{u) is not a good guess in the H ^ 0 zone, since at any finite temperature 
the local magnetizations are non zero, and therefore, messages should be non trivial even in the 
paramagnetic phase. Around the paramagnetic solution the self consistent integral equations ([8|) 
can be linearized for the first moments of the distributions q{u) and Q{U,'^U 2 ), and the critical 


temperature is defined as the moment when the system becomes singular 




In the simplest case of the Bethe approximation, only an equation like the first one in ([8|) 
appears (excluding the presence of the Q{u,ui,U 2 ) functions), and it can be solved again in terms 
of the moments around q{u}_^= S{u), or numerically as a fixed point dynamic over a population of 


fields u representing q{u) 


m. 


o 


l\l, 


-o 


u, U u 



FIG. 10. Messages in the population dynamics scheme. In order to preserve correlations, plaquette to link 
messages are updated and stored together with the parallel link to spin messages as shown in the left part 
of this figure. Each iteration step samples the surroundings of a plaquette and generates interior messages 
(see right side). 


This population sampling method can be extended to the CVM approximation in a less trivial 
way. A representation of Q{U, ui,U 2 ) with a population of fields {U, ui,U 2 ) is not enough to attempt 
a numerical solution of equation Q, since it is not trivial how to deconvolve Q in the left hand 
side of ([8]). Still, you can pretend that a random sampling of u and {U, ui,U 2 ) populations and the 
evaluation of message passing equations (j4]) and ([5]) with random local external fields hi represents a 
randomized version of the actual message passing occurring in the system, and therefore represents 
an average system. 

In the single instance case, for every update step of ([/, ui,U 2 ) messages, it is necessary to update 
simultaneously the two (ul) messages that enters in the left hand side of ([2]) 29(]. It is prudent then 
to respect this correlation in the random sampling procedure. In order to do so, a random object 
containing a {U,ui,U 2 ) and two ul is defined (see Fig. [TOl left). It is possible to generate a new 
set of {U, Ul, M 2 ) and u^ by sampling four of these objects from a large population and using them 







17 


Phase Diagram 2D Bethe 

BP, Pop Dyn, Rep-CVM 



FIG. 11. Comparing results for the Bethe approximation. In this case, the single instance results differ from 
the two average case prediction. The fact that for a given H, the critical temperature is generally lower 
than the average calculation is a consequence of the domain formation phenomena 


as the messages acting on a plaquette (see Fig. [THl right). This newly generated random object is 
added back to the population at the end of each iteration step. 

In the population dynamic method we will define the critical temperature as the moment in 
which the population of fields develop a global magnetization. This means that the balance between 
positive and negative u fields is skewed to one side, signaling the appearance of long range order. 

In Fig. [TT]we show the critical lines obtained by examine single instances of RFIM in 2D with 
the Bethe approximation (same data as in Fig. HED, and the results obtained in average case by 
both procedures: stability analysis of linearized equation and population dynamics. Remember 
that for the Bethe case these two methods attempt to solve the same integral equation. Therefore, 
had the stability analysis been performed to all orders, it would have matched population dynamics 
results. Nevertheless, for small H, a good agreement is observed in Fig. [TTJ In the limiting case 


H = 0 all results agree on the same critical temperature, which can be found accurately 


IJ]. In 


presence of disorder, on the other hand, the single instance critical line predicts for the same field 
intensity a transition temperature which is lower than the average case values. This divergence is 
closely related to the domain formation picture, because the system needs lower temperatures to 
be completely magnetized. 


The same comparison is done in Fig. [12] for the plaquette-CVM approximation. Again, at low 
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Phase Diagram 2D Kikuchi 

GBP, Pop Dyn, Rep-CVM 



FIG. 12. The GBP approximation match quite well for small H values with average case predictions. 
Galculations for the Replica CVM are valid only for weak magnetic field, whereas, for high intensities, 
approximations made in calculations break down. Error bars for the GBP graphic fall within the symbol 
shape. 


H both average case calculations coincide. Furthermore, they coincide also with the results for the 
phase diagram obtained within the single instance scenario. 


At this point it is interesting to compare these results for CVM with previous predictions for 
the EA model 3^. Unfortunately the replica stability analysis in RFIM can only be accurate in 


the low field zone, restricting our analysis to this part of the diagram. While in the 2D EA model 
the population dynamic and the replica stability give different transitions temperatures, the first 
one related to the appearance of non paramagnetic GBP solutions and the second pointing the 
lack of convergence of the message passing, in REIM both methods seems to coincide (in the low 
H region). They also coincide with the appearance of magnetized solutions in single instances, and 
asymptotically with the lack of convergence for L —?■ oo. In this sense the relation of single instance 
behavior of GBP and both average case methods is connected in a similar fashion as in 2D EA. 
The rest of the diagram (higher H) is harder to analyze, since replica stability is a priori wrong, 
and population dynamics gets into the non convergent region of single instances. 
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VI. CONCLUSIONS 

The Random field Ising model is a paradigmatic disordered model in statistical mechanics whose 
full comprehension in general dimensions resists analytic methods and poses major simulation 
difficulties. We have studied two approximations to the free energy of this model in 2D, namely 
Bethe and plaquette-CVM approximations, both in single instances and in the average case. 

The qualitative panorama of both approximations is consistent. BP and GBP on finite size 
instances and average case calculations predict transitions to an ordered phase in a region of low 
temperatures and low fields. Although it is known that this phase is thermodynamically unstable, 
the prediction of long range order in 2D is not a surprise, since mean field like approximations tend 
to stabilize ordered phases. However, Bethe and CVM differ in many aspects. The Bethe approxi¬ 
mation predicts a critical line that is above the one of the plaquette-CVM, which is consistent with 
the expectation that a more precise approximation yields results closer to reality (no critical line at 
all). On single instances the Bethe approximation converges to many solutions and GBP to only 
one or at most a few ones related by flipping large spin clusters. Moreover, our results suggest that 
GBP does not converge for large enough lattices and more important, at least for small values of 
the field ~ 0 , the non-convergence seems to coincide with the para-ferro transition temperature 
predicted by the population dynamic and the replica stability analysis. As in the 2D EA model, 
the behavior of GBP is closely connected with the average case prediction. 


VII. APPENDIX: GAUGE INVARIANCE OE GBP EQUATIONS 


GBP equations enforce the correct marginalization of the expectations (beliefs) at each level 
onto the expectation at its children levels, in a region graph approximation to the free energy of 
a model. As discussed in 3^, each marginalization requires a Lagrange multiplier, and the self 


consistency of this multipliers becomes the usual message passing algorithm. 

When implementing GBP in single instances, it makes no damage if some of the constraints are 
forced more than once. For instance, if a plaquette’s belief ^(cri, £ 72 , us, £ 14 ) is forced to marginalize 
onto two of its links b(ai,a 2 ) and 6 ((J 4 ,(Ji), then the marginals of these two beliefs must be con¬ 
sistent on variable £ 74 . Therefore, if one of them, lets say b(ai,a 2 ) is forced to marginalize on the 
belief 6 (£ 7 i), then it is unnecessary to force the other to marginalize on it, since it inevitably does. 

This means that the general prescription for implementing GBP in 36( is redundant, and the 


introduced messages or their equivalent cavity fields are under-determined. In [33| we discussed 
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FIG. 13. One of the gauge invariances of CVM approximation. Since fields are under determined, some can 
be raised(thicker arrows) and others diminished (thinner arrows) without altering the total field acting on 
a site. Eliminating this free modes by fixing some fields can help improving convergence properties of the 
algorithm. The procedure of fixing the gauge does not affect the belief or other relevant magnitudes. 


INI N I 

I I ' I 1^1 



f \ t I \ t t 

I si I si I 

I M I M I 


FIG. 14. Fields shown in bold face in this figure enter linearly into each other update equation. This fact 
can be exploited to update them in a simultaneous and coherent way. In addition to the gauge fixing this 
heuristic helped us improving convergence near the transition line. 


this fact for the Edwards Anderson 2D model, and showed a way to fix the gauge invariance in the 
cavity helds that appears as a result of their indeterminacy. In Fig. [T3]we show schematically one 
of the gauge invariant transformations of the messages. Two of the messages pointing to spin on 
bottom-left can be raised by the same arbitrary amount, and the two others reduced in the same 
amount without altering the fixed point of the equations. 

Fixing the gauge does not alter the type or the number of solutions found for the beliefs (which 
is the physically meaningful quantity, not the fields or messages). However, it can help gaining 
some more convergence in single instances, and furthermore, becomes necessary to make a correct 
population dynamics representation of the average case. We decided to hx the gauge by setting to 
zero one of the fields acting over spins in the plaquette-to-link message (C/, ui,U 2 ), as shown in Fig. 
M Furthermore, since all the bold faced fields in this figure participate in each other equations in 
a linear way, we solved the set of linear equations in each message passing step, therefore moving 


all these fields in a consistent way towards agreement 


33|]. 
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Prob Conv S.l. 2D 
H=0.5 



FIG. 15. Comparing GBP with and without gauge fixing in a sample. Probability of convergence is plotted 
as a function of temperature for a fixed external magnetic field. Although near the transition both variants 
get instable, the one with the fixed gauge converges much better. 


This mixed strategy of fixing the gauge and a consistent updating of all linearly dependent fields 
showed to improve the convergence of the GBP in the 2D Edwards Anderson model 33|. In Fig. 


we show that this is also the case in the random field Ising model. Most of the convergence 
problem appearing near the transition line from para to long range ordered phase disappears if the 
GBP is implemented with these two prescriptions. However, the convergence problems appearing 
within the long range ordered phase, did not disappear, and a growing island of non convergence 
still exists for low temperatures and fields, as discussed in the article. 
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